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Abstract 

The influence of nonmagnetic doping on the thermodynamic properties 
of two-leg S = 1/2 spin ladders is studied in this paper. It is shown that, 
for a weak interchain coupling, the problem can be mapped onto a model 
of random mass Dirac (Majorana) fermions. We investigate in detail the 
structure of the fermionic states localized at an individual mass kink (zero 
modes) in the framework of a generalized Dirac model. The low-temperature 
thermodynamic properties are dominated by these zero modes. We use the 
single-fermion density of states, known to exhibit the Dyson singularity in 
the zero-energy limit, to construct the thermodynamics of the spin ladder. 
In particular, we find that the magnetic susceptibility x diverges at T — > 
as 1/Tln^(l/T), and the specific heat behaves as C oc l/ln'^(l/T). The 
predictions on magnetic susceptibility are consistent with the most recent re- 
sults of quantum Monte Carlo simulations on doped ladders with randomly 
distributed impurities. We also calculate the average staggered magnetic sus- 
ceptibility induced in the system by such defects. 

I. INTRODUCTION 

Spin ladders have attracted considerable attention of theorists and experimentalists in 
recent years [||] . The main distinctive feature of these systems is the existence of a spin gap 



in the excitation spectrum (for ladders with even number of legs). The spin gap has exper- 
imentally been observed, for instance, in SrCu20'i systems [0. The spectrum of the spin 
ladders is rather similar to that of (integer spin) Haldane chains and spin-Peierls systems. 

Somewhat suprisingly, these gapped systems exhibit interesting behavior when doped by 
nonmagnetic impurities. In particular, Lai-xSrxCu02.5 compositions have been investigated 
experimentally and a metal-insulator transition was found 0. In the low-doping regime 
{x < 0.05) these systems show an insulating behavior, i.e., holes are localized. Later on, 
Zn doping has also been realized in the Sr{Cui^xZnx)20z compounds and a transition 
from the singlet state to antiferromagnetically (AF) ordered state was observed even at 
very low dopings ( of the order 1%) These doped systems also exhibit a substantial 
linear in T specific heat showing abundance of low-energy excitations. Very recently, the 
neutron scattering data on these doped compounds reveal a finite density of states at zero 
energies, being consistent with the specific heat data, while the amplitude of the spin gap 
itself remains almost unchanged ^ . 

Theoretically, the effects of non-magnetic doping in ladder compounds have been studied 
extensively, using various numerical techniques P-pH]], bosonization [|12|, real-space renor- 
malization group (RSRG) |]13[, nonlinear cr-model Liouville quantum mechanics ||16|| , 



the Berezinskii diagram technique [0, the supersymmetric method |T^, etc. Intuitively, 
one might expect that, because of the spin gap, the impurities would be irrelevant and 
have no significant effect. However, this is not true. The nonmagnetic impurities create 
low-energy (in-gap) localized states which dominate the low-temperature thermodynamics. 
Up to now, there is already some consensus in the theoretical understanding of this issue 
(without making explicit references) : The nonmagnetic impurity induces a spin 1/2 degree of 
freedom around it which leads to a Curie-like uniform susceptibility; the effective interaction 
between these "free" spins can be ferromagnetic or AF, depending on the impurity configu- 
rations; there should be zero-energy states which show up in neutron scattering experiments 
and give rise to additional specific heat; the AF- magnetic correlations are enhanced around 
the impurity, etc. 

However, there are still several important open questions: What is the exact form of the 
Griffiths singularity for the susceptibility (for which the RSRG analysis |13[ and quantum 
Monte Carlo (QMC) simulations [ill] g^-^e very different results)? How to derive the "zero 
modes" assumed in several calculations? How to calculate the staggered susceptibility? How 
to construct a self-consistent theory for the singlet-AF phase transition? 

In this paper we present a more systematic theoretical study of the doping effect in 
two- leg spin 1/2 ladders, based on symmetry analysis, bosonization technique and mapping 
on random mass Dirac(Majorana) fermion model. We will limit ourselves to the insulating 
regime when doped holes are localized. We first introduce the bosonization technique and 
show explicitly how the presence of holes will modify the motion of the spin degrees of 
freedom (first for a single spin 1/2 chain in Section Then in Section |T| we map the doped 
spin ladder system onto a model of Majorana fermions with random masses. Moreover, in 
Section |IV| we investigate the symmetries and the states of the fermionic model with a single 
mass kink to explicitly derive the zero modes. In Section |^ the thermodynamic functions 
are evaluated and the calculated unform susceptibility is compared with the QMC results 



11 1 . Furthermore, in Section VI we calculate the average staggered magnetic susceptibility 



caused by a defect. Finally, some concluding remarks are given in Section |V11] . 
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Some of the presented results are already known to some extent by now, especially in 
view of the similarity between the ladders and the spin-Peierls systems [|l9ip0| (although 
there are some essential differences between these two cases). However, we believe our study 
sheds new light on the problem and provides justification for several assumptions made in 
earlier papers. The comparison with QMC results shows first evidence of nontrivial physics 
in the problem-difference between " typical" and " average" behavior of the spin correlations 
in systems with randomly distributed impurities, as will be explained later in Section |V|. 



II. NONMAGNETIC IMPURITIES IN ITINERANT SPIN 1/2 CHAINS 

Let Hfjuik be a standard Hamiltonian for a one-dimensional interacting electron system. 
We shall not write Hi,uik explicitly in terms of the electron field operators (its bosonized 
version is given below), referring the reader to the literature instead We assume that 



interaction is spin-rotationally invariant, and choose the interaction constants in such a way 
that the spin sector of the model remains gapless, while charge excitations have a finite gap, 
a repulsive half-filled Hubbard model being a typical example. Although we will only be 
dealing with the spin degrees of freedom in this paper, it is very helpful to include explicitly 
the charge degrees of freedom as will be clear from the later presentation. Then, at energies 
well below the charge gap, Hi,uik describes an itinerant SU(2)-symmetric spin-1/2 chain. 
The electron field operator is bosonized as 



'RJx)+e-'''^''LJx) 



with 



RJx) 



LJx) 



Oq being the short-distance cutoff. The chiral Bose fields (pR^L)^ combine in the standard way 
to produce the charge phase field $c and the spin phase field as well as the corresponding 
dual fields Qc{s)- As usual, at low energies, the charge and spin degrees of freedom decouple in 
the bulk: Hh^ik = Hc + Hg. The charge Hamiltonian has the form of a quantum sine-Gordon 
model 



vrao 



dx cos 



(1) 



where rric > is the bare charge mass, and the phase field is rescaled according to $c 
\/TQ^c, with Kc < I being the charge exponent. Here Hq is the canonical Hamiltonian for 
the Gaussian model 



dx 



where He is the momentum canonically conjugate to $c, and Vc is the (charge) velocity. Up 
to a marginally irrelevant perturbation, the Hamiltonian for the spin degrees of freedom is 
simply 
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(2) 



with the appropriate spin velocity Vg. 

In the rest of this Section we discuss the effect of impurities in a single chain. A single 
nonmagnetic impurity in spin-1/2 chains has been studied, for localized spins, by Eggert 
and Affleck (EA) [^. Let us start by summarizing their results. EA discovered that there 
are impurities of two types: 

(L) Impurities which violate the site parity P5. These impurities may or may not respect 
the link symmetry P^. An example is an exchange constant modified on a single link. Such 
impurities are relevant and break the chain up. The infrared stable fixed point corresponds 
to an open boundary condition. 

(S) Impurities which respect Ps (hence violate Pl)- An example is two neighbouring 
exchange couplings modified by the same amount. These impurities are irrelevant, so that 
at low energiers the chain "heals" . 

The physical interpretation of EA's findings is as follows. Given that the spin dimer- 
ization is the leading instability of the Heisenberg spin-1/2 chain, one can immediately see 
that the dimerization order parameter can be locally pinned by the L-type impurities but 
not by the S-type impurities. Another way around is to say that the dimerization operator 
e(x) is invariant under P^ but changes its sign under P^: a local relevant perturbation e(0) 
is allowed for L-type impurities, but it is, by symmetry, prohibited for S-type impurities, 
with 5a;e(0) being the leading irrelevant operator. 

In addition to EA's considerations we need to trace the spatial behavior of the charge 
phase field ^d^) in order to understand how the coupling between the chains is modified by 
the presence of these irrelevant impurities. (We imply S-type impurities: see the discussion 
at the end of the Section.) 

Let us first consider the case when the charge sector is gapless. The system admits an 
arbitrary electron charge induced by the impurity potential: 




For an impurity localized at the origin over a scale ~ it means that 

9{x) being the step function defined by 9{x < 0) = 0, 6'(0) = 1/2, and 9{x > 0) = 1 [in the 
gapless case, an arbitrary constant can be added to Eq.(D]. Eqs. (D and (D describe the 
well-known charge fractionization in a Luttinger liquid p3| . 

On the other hand, when the charge sector is gapful (the case we are really interested 
in), the ground state expectation value of the charge phase field should coincide with one of 
the Zoo-degenerate vacua of the potential in the sine-Gordon model (|l|): 



($c(x)) y ^ ^ ' (^) 

n being an integer. It is important to notice that no local impurity potential can overcome 
the bulk energy of the system. Hence the asymptotic condition (^) must be satisfied for 
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any impurity scattering operator. The integer n, however, can vary and may be different 
for X — oo and x +00; this only involves a local alternation of the umklapp scattering 
term. 

The condition has an immediate effect on the value of the electron charge @ that 
can be trapped by impurities in a gapped system. Indeed, comparing (^) and (||) one finds 
that the charge is quantized 



Q = em, 



(6) 



with m being another integer number equal to the increase of n in going from x = — 00 to 
X = +00. 

For a single impurity it is natural to assume that m = 1 (m = — 1), so that one electron 



(hole) is trapped around the impurity site ||2^. If many such impurities are scattered along 
the system, then the charge phase field acquires an average value 



2K, 



■N{x) , N{x) = Y,0{x- 



(7) 



The bosonized form of the staggered magnetization is then given by [p5 



Tfix) 



X{x) . 
sm 



Trao 



n^ix) 



exp 



with the function \{x) defined as 



\{x) = (cos 



27T Kr(!>r(x) 



Using (|^), one obtains 

A(a;) = Aoexp [inN^x)] = Ao]^sign( 



(8) 



(9) 



(10) 



where Aq is a non-universal dimensionless constant equal to the average d) in the absence 
of the disorder. 

Similarly, the spin dimerization operator acquires an identical x-dependent prefactor 
due to the alternation of the average value of the charge phase field by the nonmagnetic 
impurities: 



e{x) = (-1)"S„ ■ S„+i 



A(x) 



■ COS 



^<i>s{x) 



(11) 



Thus, the only but important effect of S-type impurities on the single chain is the appearance 
of the sign alternating factor in the definitions of the staggered component of the spin density 
and dimerization field. The consequences of this factor for the model of coupled chains are 
discussed in the next Section. 

Here we would like to note that Sr doping of Ref. ^ probably leads to impurities of 
the type S. Indeed, Sr doping produces holes in the system. Given the similarity of the 
chemical composition of La — Cu — O chain systems and LaCu204^ high-T^ compounds, it 
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is natural to assume that holes are localized on the oxygen ions (for they are known to be 
localized on the oxygen ions in the high-T^ materials). The localized hole with spin 1/2 
represents then a "new" site in the magnetic chain. Since this effectively adds an extra site 
as compared to the pure chain, the physical meaning of the staggered magnetization sign 
change becomes almost trivial. Since S-type impurities are irrelevant in the renormalization 
group sense (i.e., the coupling of the hole spin to its neighbours flows toward the uniform 
exchange J), this change of the sign is the only effect of Sr doping. 

On the contrary, the impurities of L type (like Zn doping of Ref. 0) are relevant, so that 
one may conclude that the chain segment model must be used. We would like to point out 
that the last conclusion is not always correct. It is true that L-type impurities are relevant 
with respect to the single chain ground state, but they are still irrelevant with respect to 
the ladder ground state since the latter has a gap. Thus, for low concentrations, the L- 
type impurities should have the same effect as S-type impurities, the crossover to the chain 
segment model occuring only at higher concentrations. (The crossover concentration is, of 
course, exponentially small in the ratio of J to the spin gap.) Therefore, generically speaking, 
one should not consider severed chains as a realistic model for nonmagnetic dopings in spin 
ladders. 



III. DIRAC AND MAJORANA FERMIONS WITH RANDOM MASS 

Let us now consider a model of two weakly coupled S = 1/2 Heisenberg chains - a two-leg 
spin ladder. Its Hamiltonian, 

H = Ho + H^, (12) 
consists of two terms. The first term, 

i=i,2 

describes two decoupled chains (the subscript s for the 'spin' is suppressed since we shall 
only work with the spin phase fields in what follows). The second term 

H± = J dx[J±ni{x)n2{x) + Uei{x)e2ix)] , (13) 

is responsible for the interchain interaction. Note that only the relevant interactions are re- 
tained while all the marginal terms, leading to renormalization of the masses and velocities, 
are neglected. J± is the interchain exchange coupling. The second term in (p!3D , which cou- 
ples dimerization order parameters of different chains, can either be effectively mediated by 
spin-phonon interaction |2^ or, in the doped phase, generated by the conventional Coulomb 



repulsion between the holes moving in the spin correlated background ||2] 

Using the bosonization formulas (P) and ( pTD and introducing the symmetric and anti- 
symmetric combinations of the spin fields. 
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one easily finds that, like in the case of a pure spin ladder p8[, the total Hamiltonian ([I2D 
factorizes into two commuting parts 



H = H. 



with H± given by 



u 



vrao 



dx mix) cos 



47r<l>+(a;) 



(14) 



and 



u + l 



-mix) cos 



2m(x) 

H ^cos 

vrao 



ATrQJx] 



47r<l>„(a;) 

}■ 



(15) 



Here u = U/ J± and 

m{x) = mt{x), m 



27T 



t{x) 



exp 



x) 



i=i,2 



(16) 



The function t{x) = ±1 changes its sign whenever one passes through a position of an 
impurity on either chain. If one assumes, as we shall do, that the positions of the impurity 
centers are uncorrelated, then the function t{x) represents a random 'telegraph process' 
characterized by the average density no of the impurities, a = l/riQ being the mean distance 
between them. Thus the disorder manifests itself in multiplying the interaction term by the 
telegraph signal factor t{x). This leads to a random mass fermion problem, as we shall see 
shortly. 

Since the scaling dimension of the interaction terms in the Hamiltonians H± is equal to 



1, these can conveniently be re-fermionized as in the pure case p5 



Let us start with H+. The chiral components of $+ can be related to chiral component 
of a new Fermi field ip by the standard bosonization formula 



exp 



±iV47r0+,i?(L)(a;) = V2vrao tpR{L)ix) . 



(17) 



In terms of this new Fermi field, becomes 



dx\^-Ws V^|j(a;)9x.V'j^(x) - ^i(x)9a;V^^(x) 



imtt{x) [V^|j(x)^i(a;) - V^i(a;)^^(x)]] } 



(18) 



where = (1 — u)m. Thus, the Hamiltonian describes Dirac fermions with a ran- 
dom (telegraph signal) mass. It is sometimes convenient to separate the real (Majorana) 
components of the Fermi field operator 



XI 



1 

V2 
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so that 



Hom= E ^m[C1 (19) 

a=l,2 

with 

^m[C] = - / [CR{x)d,U^) - Cl(x)9.Cl(x)] 

+ ^mt(x)eii(x)a(a;)} (20) 

standing for the Hamiltonian of the random mass Majorana field. 

The Hamiltonian H_ admits a similar re-fermionization procedure based, analogously to 
(p!7|), on the introduction of the Fermi field 



exp 



The only difference with (|T8|) is the appearance of a 'Cooper-type' mass due to the presence 
of the dual field (0-) in the interacting part of Eq(|l^). This can easily be diagonalized 
by directly passing to the Majorana components 

Xr{l){x) = [Cl(L)ix) + iCR(L)ix) ■ 

As a result the Hamitonian H_ decomposes into two random mass Majorana models 

H- = HT/[C']+HTAC'] 

with different amplitudes of the mass telegraph signal, rrit and = — (3 + u)m. 
The total Hamiltonian can now be represented in the form 

H = H^ + H. = j2HZ^[Q (21) 

a=0 

{niQ = rus, mi^2,3 = ^t) which, as in the pure case [^,^, reflects the spin rotational 
symmetry of the problem (SU(2) symmetry is, of course, preserved by the disorder). The 
Majorana fields C,"" with a = 1, 2, 3 correspond to triplet magnetic excitations, while 
describes a singlet excitation. All these fields acquire random masses due to the presence 
of disorder. (One should bear in mind that, when marginal interactions are included, the 
masses rrit and get renormalized, and the velocities in the triplet and singlet sectors 
become different: Vg ^ Vt, Vg.) 

The mass kinks create low-energy states within the spin gap. One therefore expects 
the low-temperature thermodynamic functions to be dominated by these states. Before 
constructing the thermodynamics of the system (Section 0), we investigate in detail the 
theory with an isolated kink. 
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IV. ZERO MODES AND FERMION NUMBER IN A GENERALIZED DIRAC 

MODEL 

Let us consider a Dirac Hamiltonian H with the structure of iJ_ in Eq. assuming 
that both the Dirac (CDW-hke) and Cooper mass functions, mi(x) and m2(x), are of the 
"telegraph process" type but otherwise are independent. It is convenient to make a chiral 
rotation 

under which only the kinetic energy term in (|15D is modified (75 

H{x) = -w{R^d^L + L^d^R) + imi{x){R^L - L^R) 

+ im2{x){R^L^ - LR). (23) 

We know that this generalized Dirac Hamiltonian factorizes into two decoupled massive 
Majorana fields: 

HmA^) = -tvCRd^a + ^m^ix)eReL, = ±)- (24) 

where 

m±{x) = mi{x) ± m2{x). 

Using this correspondence, we would like to understand under what conditions the zero 
modes of the original Dirac model p3| ) can appear, and what are the corresponding fermion 



numbers fSQ] 



Let us first make a few remarks on the symmetry properties of the model (^) and some 
of their implications. For m2 = 0, mi 7^ 0, the Hamiltonian is invariant with respect to 
global phase transformations of the fermion fields, R — > e*"i?, L —>■ e*°L, which amounts 
to conservation of the total particle number 

N = dx tp"^ {x)tjj{x). 



The continuous chiral symmetry R e'^'^R, L e~'^'^L, is broken by the Dirac mass, and 
the current 

J = dx ijj'' {x)'y5ip{x) 



is not conserved. Since N is conserved, the existence of zero modes for a solitonic shape of 



mi (a;) will lead to the appearance of fractional fermion number ||3CI|| . 

On the other hand, for mi = 0, m2 7^ 0, the global phase invariance is broken by the 
Cooper-mass term, and N is not conserved. So, there is no nonzero fermion number in 
this case. However, there is a continuous chiral (or 75) symmetry which is respected by the 
Cooper mass term, giving rise to conservation of the current J. Again, if zero modes exist, 
a fractional local current will appear. 
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This can be understood in two ways. The Cooper mass can be transformed to a Dirac 
mass by a particle-hole transformation of one chiral component of the Dirac field, e.g. R 
R, L ^ L'^ . Under this transformation, N ^ J; hence the quantized fractional current. The 
other explanation is physical. For a solitonic shape of the Cooper mass m2{x), the existence 
of a nonzero local vacuum current is a manifestation of the Josephson effect, because when 
passing the impurity point the gap function changes its phase (by vr). 

Consider now a more general situation, relevant to our discussion of disordered spin 
ladders, when both mass functions mi{x) and m2{x) have single-soliton profiles with co- 
inciding centers of the kinks, say, at x = 0, but with arbitrary signs of the corresponding 
topological charges and arbitrary amplitudes at x = ±00. Before starting calculations, we 
can anticipate the characteristic features of the solution. At mi(x),m2(x) 7^ none of the 
above mentioned continuous symmetries survives, i.e. neither N nor J is conserved. There 
are only discrete symmetries present: one is the E ^ —E symmetry generated by transfor- 
mations R ^ R, L — i> —L, the other being charge conjugation or particle-hole symmetry: 
R ^ R\ L L'^ : H H. These symmetries imply that, if localized vacuum states 
appear in the solitonic backgrounds of the masses mi{x) and 777,2(0;), those are supposed to 
be zero modes, i.e. localized states exactly a.t E = 0. Moreover, for massive Dirac fermions 
(7772 = 0), it is the charge conjugation symmetry that implies quantization of fermion number 



(see review article ||30[): (A^) = n/2, n & Z. For our Hamiltonian (23) the first statement 
still holds: if localized modes exist, they should be zero-energy modes. But the second 
statement is no longer correct: since and J are not conserved, their local expectation 
values will not be quantized. Instead, the magnitude of the fermion number or current will 
depend on the ratio of the amplitudes m\/m^ of the corresponding mass functions. Only in 
the limiting cases m\/m^ — > 00 and m\/m2 — >■ will the universal quantized values of the 
charge and current be recovered. 

Let us now turn to calculations. We consider a single massive Majorana field described 
by the Hamiltonian ( |2^ ) for a fixed index a. It can be represented as 

Hm = ^em, (25) 



where 



Ti. = pvai — m{x)(T2 

_ ( 0, pv + im{x) 

\pv — im{x), 



(26) 



Here di and (T2 are the Pauli matrices. Introducing a Majorana 2-spinor 

Ur{x) 

ul{x) 



u{x) = ( ^^f^) ] (27) 



we get a pair of first-order equations: 



{v- — \- m{x))uji = iEuL, {v- m{x))uL = iEur. (28) 

dx ax 



Let ?77(x) have a step-like jump at a; = 0: 
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m{x) = mosign x. (29) 

For niQ > we shall call configuration (pOD a soliton (s); the case mo < corresponds to an 
antisoliton (s). It immediately follows from Eqs. (^) that in the case of a soliton 

Us{x)= ( J ) uo{x), (30) 



whereas in the case of an antisoliton 



u-^{x) = J j "o(a;). (31) 

Here 

uo(a;) = Ao^^^exp(-|a;|/Ao), Aq = f/mo, (32) 

is the normalized zero-mode wave function for a bound state of the Majorana fermion, 
appearing at the discontinuity point of m(x). 

Now we have to single out the contribution of zero modes in the spectral expansion of 
the Majorana field operator: 

U^)=d[]^^u,{x)+Ux), (33) 
Ux)=d[\'^Uo{x)+h{x) (34) 

Here d is the Majorana operator for the zero mode, and ^ (x) is a contribution of the contin- 
uum of scattering states which, due to the E —E symmetry, do not affect the expectation 
values of fermion number and current. This part of the Majorana field operator will not be 
considered below. 

Let us now turn back to the Hamiltonian (0). The fermion number and current operators 
can be expressed in terms of the Majorana fields as follows: 

N=tia^R+a^L)=^iefr, (35) 

4 = + ^Ur) = Vr, (36) 

(remember that 75 = ai). There are two qualitatively different cases. 
1) Both m+ and m_ have solitonic (or antisolitonic) shapes. 
Choose, e.g. m']_,m° > 0. Then 

^+(x) = d+ 1^ J j M+(x), 

= d_ 1^ jyuo(^)- 

In this case one finds that the local current vanishes, while the zero-mode fermion number 
is given by 



11 



N= {N)= I dx (: ^Ij\x)^Ij{x) :) = 7-^d+d_ 

1 1 



cosh 9 

Here we have introduced a local complex fermion operator 



a'a — 



(37) 



V2 ' 

while 

exp 26 = m°/m°, (38) 

(This parametrization assumes that both m*^ and m° are positive, i.e. > JTig.) The factor 
1/ coshB represents the overlap integral between two zero-mode wave functions in the (+) 
and (— ) channels, respectively. One finds that for a pure Dirac mass = 0, = m° 
and G = 0. As a result, the fermion number is quantized: 

iV = ata-i, (iV) = 4- 

2) m_|_ has a solitonic shape, while m_ has an antisolitonic shape, or vice versa. 
We shall choose m° > 0, mP_ < 0. Then 

^+{x) = d+ 1^ J j <(x). 

In this case the situation is inverted: the vacuum fermion number vanishes, while the current 
does not. The same calculations lead to 

J = —^ia^a--). (39) 
coshe^ 2' ^ ' 



Let us now turn to the Hamiltonian if_ in Eq. (p^. In this Hamiltonian m+(a;) = 
3m(x), m_(x) = — m(x), so that, for m[x) = mo sign x with arbitrary sign of mo, we are 
dealing with the case of a nonzero local current (case 2)). Since cosh 6 = 2/^3 m our case, 
we get: 

J=^(ata-1). (40) 

These simple results can be immediately applied to the smooth parts of the spin density 
at the impurity point. One easily finds that the total spin density 

SUx) = SJ(x) + S^{x) = -d^<i>+{x) =: ij\x)ij{x) : (41) 

TT 

while the realtive spin density 
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Sl{x) = Stix) - = -d,^-{^) =: x\x)x{x) : (42) 

So, {S^{x)) and {S^{x)) coincide with fermion numbers in the (+) and (— ) channels: 

{SUx)) = ^-crul^{x), a = ±l), (43) 
{Sl{x)) = 0. (44) 

Notice that the fermion number is quantized; hence the locahzed spin is 1/2 (as can be 
checked by integrating (S'^(x)) over x). 

The spin current operators are defined as 

]{x) = v[Jn-3L], (45) 

so that we have: 

{jl{x)) =V a U-rno{x)u3mo{x), (46) 



hence the integrated current 



0-) = ^va. (47) 



Therefore the spin current at the impurity point is not universal and it is determined by the 
ratio of the singlet and triplet masses. 



V. THERMODYNAMIC FUNCTIONS FOR DISORDERED SPIN LADDERS 

A. Free energy 

Since we shall be interested in the behavior of the system in a magnetic field, we start 
this section by adding the magnetic field term to the Hamiltonian (in this paper only a 
spatially homogeneous magnetic field will be considered): 

-fisH J dx [slix) + s',{x)] = I dxd,<!>+{x) 



As follows from ([T7|) , the magnetic field appears as a chemical potential (equal to ^bH) in 
the refermionized version (PI^D of the if+ part of the Hamiltonian 



- ^ibH j dx : \'4)^R{x)i)R{x)+i)l{x)i)L{x)\ ■ (48) 



(the normal ordering in this formula must be taken with respect to the ground state of the 
system without magnetic field). 
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The magnetic field breaks the spin rotational invariance of the problem. Hence it is 
convenient to work with the Dirac version (^) rather than the Majorana version (|T9|),(^0D 
of the Hamiltonian We thus rewrite the total Hamiltonian for the spin ladder in the 
magnetic field as 



H = H 



mt;H ^ 



+ HZ/[C'] + K/[C' 



(49) 



Assume that we know the exact, averaged over the disorder, free energy for the random 
mass Dirac fermions {H^'^[ip]) as a function of the magnetic field H and temperature T. 
This free energy, which we denote by FJ^{T, H), will be actually calculated in what follows. 
Clearly the free energy of random mass Majorana fermions F^{T) can be inferred from the 
above function (notice that there is no magnetic field term for a single Majorana fermion). 
Since H^'^~^ decomposes into two Hamiltonians H^^ for independent Majorana fields and 
the free energy is an extensive variable, we simply have 



F^{T) = -F^{T,0) 



The total free energy of the system (^9|) is therefore 



F(T, H) = F™' (T, H) + ^F^^ (T, 0) + ^F™^ (T, 0) . 



(50) 



In the following we shall focus on the function FJ^{T, H). Since the fermions are nonin- 
teracting, the regularized free energy can be written as 



AF™(T, H) = F^{T, H) - F™(0, 0) 





(51) 



T / depD(e)ln 



1 + e' 



T / depD(e)ln 



1 + e' 



where (3 is the inverse temperature, and Pd{^) is the density of states for the Dirac fermions 
(|I^) , averaged over the disorder |^ . 

In fact, PD(e) is a single particle density of states for the quantum mechanical problem 
of a Dirac particle with a random mass [^. The wave-function of the latter, for which we 
keep the same notation as for the field operator 



satisfies the Dirac equation of the form 

[-ivtasdx + m{x)(72] = eip 



(52) 



[there is no chemical potential term in Eq. (p2D, for it has been explicitly separated in (pl|)]. 
It is convenient to make the chiral rotation 



V2 



[v{x) ± u{x)] , 
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The spinor component u{x) then satisfies the Schrodinger-type equation 



mi 



x] 



u 



Eu 



(53) 



where fh{x) = m{x)/vt and E = e'^/vf. The spinor component v{x) satisfies Eq. (53) with 
m{x) replaced by — m(x). 

The equation (^) is known as Witten's toy model for supersymmetric quantum mechan- 
ics |^B|]. To our great advantage, this problem with a telegraph signal m(x) has recently 
been solved exactly by Comtet, Desbois, and Monthus (CDM) In particular, CDM 

calculated the disorder-averaged integrated density of states N{E) for the problem (^3|). 
Comparing Eq.(|53|) and Eq. ([5^) , one easily finds that the density of states for the Dirac 
particle we are interested in is related to the integrated density of states of Ref. |3^ by 



PD[e) 



2e 



N' 



(54) 



Thus we have, in principle, an exact solution for the free energy of the disordered spin 
ladder. The analytical expression for the function N{E) found by CDM is, however, quite 
complicated. Therefore, instead of reproducing this expression here (an interested reader is 
referred to the original paper p5), we shall consider particular limiting cases. 



B. Low energy thermodynamics 

In this section we consider the thermodynamic functions of the disordered spin ladder 
at the lowest energy scale, e -C eo, with eo to be determined later. 

According to the CDM solution and Eq.f^lf), in the e — > limit, the density of states 
takes the form (see also discussion in the next section) 



where 



o-Q = 4^ • (56) 

vino 



The expression (^5]) is given in the leading logarithmic approximation. 

Using (pT]) together with the particle-hole symmetry of the problem, implying that 
Pd{^) = Pd{—^), the magnetic moment of the system, M(T,H), induced by the external 
magnetic field is found to be 

oo 

M(T, H)=fiBj depoie) [f{e - /i^^) - /(e + /i^//)] , (57) 



where /(e) = e^*^ + 1 is the Fermi function. 
The linear magnetic susceptibility therefore is 
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XliT) = ^ J depoie) cosh^ 



-■^J. (58) 



As the temperature is lowered, the hnear susceptibihty diverges as 

This can be interpreted as a Curie hke behavior with a vanishing Curie constant, C{T) ~ 

On the other hand, the zero-temperature magnetic moment is simply proportional to 
the integrated density of states 



For a weak magnetic field 



(60) 



leading to a nonlinear susceptibility 



(^) = ^(0>^) ^ /^^^o (62) 



which diverges in the same way as the linear susceptibility (|59| ) does. This indicates that 
the magnetic field scales as the temperature, as it should be for noninteracting particles. 
The differential susceptibility, however, directly measures the density of states 

X.{H) = = f^BPoif^sH) ^ ^^^3(,/^) • (63) 

Interestingly, the low-temperature correction to the magnetic moment in a finite field is also 
a highly singular function of the field 

MiT.H)-MiO.H).--j;f^^. (64) 

The zero-field free energy reads 

oo 

AFd{T, 0) = -2T J depDie) In (l + e"'^^) . (65) 



The low-temperature free energy correction is therefore 

AF,(T,0)^-^^^. (66) 
^ ' ' In2(l/T) ^ ^ 

An unusual behavior of the specific heat follows: 
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(67) 



where the parameter 



V 



is a function of the ratio of the velocities in the singlet and the triplet sectors. 
The low-temperature entropy vanishes as 

^^^^ " ' 

indicating a non-degenerate ground state. However, the entropy vanishes very slowly with 
temperature reflecting the presence of 'quasi-free moments' in the system (see also the next 
Section). 

According to (|67|), the specific heat coefficient, Cv(T)/T, diverges as T — > 0. Yet the 
specific heat coefficient is less divergent than the linear magnetic susceptibility (^). Hence 
a very large Wilson ratio 

We notice that if one instead associates the Wilson ratio with the differential magnetic 
susceptibility (§5), then this modified Wilson ratio becomes a constant, only depending on 
the ratio of the velocities: 

^^^^^ = Cy{T) = h[23^ 



C. Intermediate regime 



The singularity in the density of states of the form (|55|) has been obtained by Dyson 
back in 1953 |^ for a model of disordered harmonic chain. In the electronic spectrum at 
the centre of the Brillouin zone such a singularity was identified by Weissman and Cohan 
| 36| for the case of a non-diagonal disorder (random hopping model). The latter model is in 
fact directly related to the random mass Dirac problem through the notion of zero-modes 



(see Section |I^ and also below). The Dyson singularity in the density of states persists the 
case of a half-filled electron band with random backscattering, as shown by A. A. Gogolin 



and Mel'nikov [|3^ who also obtained the low-temperature asymptotics for the magnetic 
susceptibility (|59D to explain experimental data on TCNQ salts. A similar low-temperature 



magnetic susceptibility has been predicted by Fabrizio and Melin for inorganic spin- 
Peierls compounds CuGeOs. It must be pointed out that the spin-Peierls systems, sharing 
with the spin ladders the property of having the spin gap, behave in quite a similar way 
under doping 
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The low energy behavior (|55D is characteristic for various particle-hole symmetric models 



of disorder and most probably represents a universality class. For the random mass Dirac 



problem such a behavior was found by Ovchinnikov and Erikhman assuming a Gaussian 
white noise distribution of the mass variable m{x). The advantage of the CDM solution 
for a telegraph signal mass, which incorporates the Gaussian distribution as a particular 
case, is that it keeps track of the high-energy properties extrapolating between the universal 
low-energy regime (e -C eo) effectively described by the Gaussian distribution, and the high- 
energy regime (e ^ mo) of almost free massless particles. This enables us to describe the 
intermediate regime (eo ^ e < mo). The latter only exists for low impurity concentration. 
Indeed, CDM found that for no — > the integrated density of states, after an initial increase 
at low energies, saturates to the value 

A^(e) ~ y at eo < e < mo . (71) 

Consequently, the density of states Pd{^) almost vanishes in this region. From ([55| ) and (|7TD 
we can roughly estimate the crossover energy as 

eo ~ moexp ( ^ ^^2) 
V vno J 

Let us now consider the disordered spin ladder at temperatures eo <C T < m. From 
Eq.(^8l) for the magnetic susceptibility we obtain 

MT) ^ ^ . (73) 

This is exactly equal to a magnetic susceptibility of free spins S = 1/2 with concentration 
uq. This is in agreement with our discussion of zero modes in Section Eqs (|59|) and (|73D 
mean that the Curie constant, being almost temperature independent for eo <^ T < mo, is 
quenched in the region of temperatures smaller than eo. This behavior is different from the 



one found by Sigrist and Furusaki [|I^ who, in particular, predicted a finite Curie constant 
at T ^ 0. We will compare in Section [VD| the recent QMC results [|lT| with theoretical 
predictions which seems to show that our prediction is confirmed. 

It is instructive to consider the free energy correction (35) at eo ^ T < m, from which 
it follows that the entropy 

^(r)~21n2no. (74) 

Bearing in mind that a local Majorana fermion has a residual entropy In ^/2, we conclude that 
the expression (|74|) indicates the presence of four local Majorana fermions at each impurity 
location. Clearly three of these local Majorana fermions originate from the bulk triplet 
mode, while the remaining one is due to the singlet mode. An effective Hamiltonian for the 
local Majorana fermions (which we denote as ^* ) can be written on purely phenomenological 
grounds. Indeed, the effective Hamiltonian, that respects all symmetries of (^^, takes into 
account the fact that the magnetic field couples to the a = 1, 2 components of the triplet 
mode, and finally preserves the quadratic nature of the problem, must be of the form 
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H, 



E 



E 



-t- T^jUqUq 



HBH^d\d: 



2 ) 



(75) 



where r*^*^ are the random hopping integrals related to the overlaps of the zero modes in the 
triplet (singlet) sector, while c?* are the individual Majorana zero-modes operators studied 
in Section The expression (|75| ) clarifies the relation of our random mass problem with 
the non-diagonal disorder problem of Ref. 



36] and hence with the original Dyson model 



|35|. 



D. Comparison of theory with quantum Monte Carlo simulations 

We have derived above formulas for various thermodynamic quantities at low tempera- 
tures. It would be very important to check these formulas directly by experiments. Unfor- 
tunately, these experiments are very difficult to perform: one reason being the interference 
of various factors (lattice, other impurities, etc); another reason being the smallness of the 
logarithmic corrections. On the other hand, there has been just performed a set of very nice 



QMC simulations on doping effect in two- leg ladders [^. Those authors could simulate 
very large systems (up to 2000 sites) and get down to very low temperatures (T = 0.005 J). 
In particular, they could carry out random average over different impurity configurations 
(up to 20 realizations) which is a crucial factor in comparison with analytic theory (which, 
of course, assumes random distribution). 

In Fig.l we replotted their numerical data on uniform magnetic susceptibility (Fig. 2 
in their original paper) vs l/(lnT)^ (The temperature T is measured in units of J). The 
fitting formula is xT = c [ a + 6/(lnT))^], where c is the doping concentration, while 
the parameters a = 0.185(3), 0.145(2), 0.126(1), and b = 0.43(3), 0.29(1), 0.23(1) for c = 
0.01, 0.05, 0.1, respectively. It is important to note that according to RSRG considerations 
]r3| a should be 1/12 at very low temperatures, and 1/4 at intermediate ones (as indicated 
by arrows and dotted lines in Fig.l), whereas b should be zero. On the contrary, according 
to our analysis, eq. (|59|) a should be zero, while b should be 1/2. The numerical results 
do clearly show the presence of the logarithmic term, but the effective Curie constant does 
not vanish entirely as T — > 0: There is a finite intersection a ^ 0, and the slope b is less 
than 1/2 as expected. As mentioned in the previous subsection, the asymptotic logarithmic 
behavior should be valid for T < eo, where eo is the low energy scale in the problem (roughly 
speaking, the soliton band width). Therefore, one should anticipate good agreement only 
at very low temperatures ( when l/(lnT)^ <^ 1). It is quite interesting to notice that the 
linear fitting is better (over broader temperature range) for higher concentrations, where eg 
is larger. 

Our tentative interpretation for the absence of full agreement with the theoretical pre- 



diction is due to the fact that the random sampling in |11| is still not big enough to fully 
demonstrate the anticipated behavior. To explain this point, let us recall the basic physical 
picture of the zero energy states in the Dirac model with random mass. Some of these 
states are genuinely localized, while the others are only "quasilocalized" . The first category 
of states is "typical", while in taking random average the second category of states does 
contribute in a substantial way. These features show up clearly in the spin-spin correlations 
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calculated using the Liouville quantum mechanics [|T6[, the Berezinskii diagram technique 



T7| , and the supersymmetric method |]T8[. The typical configurations for the spin correlation 
functions are exponentially decaying, whereas the average behavior has a power-law decay. 
The difference between the "typical" and "average" configurations is the nontrivial piece 
of physics involved in randomly doped gapped spin systems (as well as some other random 
systems). The fact that the density of states shows the Dyson singularity and other thermo- 
dynamic quantities show logarithmic singularities are all due to the same reason. It is quite 
remarkable that the signature of this behavior has shown up in the QMC simulations. As 
clear from the above explanation, the logarithmic singularity will show up as " full-fiedged" 
only if the random sampling is really big. Otherwise, we will still see some constant term 
a as "remnant" feature of the dominance of exponentially decaying states. Hopefully, with 
the further improvement of the numerical techniques, this prediction can be checked more 
precisely. Namely, when the sampling becomes bigger and bigger, the intersection with the 
vertical axis (the remaining part of the Curie constant) should gradually vanish. To the 
best of our knowledge, there is no explanation for this type of logarithmic singularities other 
than the one described above. Therefore the presence of this term per se in the numerical 
simulations is already significant. 



VI. STAGGERED MAGNETIZATION NEAR THE DEFECTS 

In fact, in the continuous Majorana model there are two vacuum averages: the staggered 
magnetization and the smooth magnetization are both nonzero in the vicinity of the point 
where m{x) changes sign. Unfortunately, we have not been able to calculate the staggered 
magnetization for the model with sign-changing m(x); we have done it only in the model 
with a sharp edge (that at the end of a broken chain). Nevertheless since this solution shows 
the presence of zero modes we think that it gives a qualitatively correct description of the 
staggered magnetization. 

The calculation is based on two facts. The first one is that the order and disorder 
parameters in the Ising model are expressed in terms of fermion creation and annihilation 
operators R{e), R+{e) as follows (T > |g: 



/i(r,x) =: exp[^pir(r,x)] : a{T,x) =: ?/;o(r, x) exp[^pi.(r, x)] : 

/oo 
d9id92imh[{ei - 92) /2] exp[(^i + 92) /2\ 
-00 

X exp[— mr (cosh 6'i + cosh 6*2) — imx{smh.9i + sinh ^2)]-R(^i)-R(^2) 

-|-terms with (76) 

/oo 
d9{e^''^ exp[— mr cosh 9 — imx sinh 9]R{9) 
-00 

-hterm with (77) 
These fermion operators satisfy the standard anticommutation relations: 
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[R{e),R-^{6')]^ = 5{6-6') 

which, in the case of the Ising model represent the simplest realization of the Zamolodchikov- 
Faddeev algebra. 

Since the operators are normally ordered and < 0|-R~'"(^) = in < 0|/i and < 0|(j we can 
omit all terms with R^: 

< 0\fi{T,x) =< 0|exp{— / (i6'id6'2tanh(6'i2/2) exp[(6'i + 6'2)/2] exp[-mr(cosh6'i+cosh6'2) 

2 J—oo 



- imx{smhe^ + sinh^2)]i?(^i)i?(^2)} (78) 

/oo 
(i6'[e^' ^ exp[— mr cosh 6 — imx sinh 6]R{6) 
-oo 

— i 

X exp{ — / c?6'i(i6'2tanh(6'i2/2)exp[(6'i + 6'2)/2] exp[-mr(cosh 6'i + cosh6l2) 

2 J—oo 

-zma;(sinh^i + sinh^2)]i?(^i)i?(^2)} (79) 



The second fact is that in the approach suggested by Ghoshal and Zamolodchikov ||41 
time and space coordinates are interchanged and the boundary is thought about as an 
asymptotic state at t — oo. This out-state is denoted as \B >. Each integrable model has 
its own \B >- vector. For the Ising model with free boundary conditions can be represented 
by the state vector is given by 

? /"OO 

\B >= [1 +i?+(0)]exp{-- / decoth{e/2)R-^{-e)R^{9)}\0 > (80) 

2 J — oo 

Notice that it contains a fermionic creation operator with zero rapidity; this operator cor- 
responds to the Majorana zero mode - a boundary bound state. This mode is the zero 
energy described in the previous sections. Since in the Ghoshal-Zamolodchikov formalism 
space and time are interchanged, -R^(O) formally creates a state with zero momentum. 
Let us calculate a vacuum average of fi at point X: 

< fi{t,X) >=< fi{T = X,x = t)\B > 

/OO j poo 

d9ide2A{ei,92)Ri9i)Rie2)]exp[ — / d9 coth{e/ 2) R+{- 9) R-^ (9)} \0 > 
-oo 2 J—oo 



oo 

oo 



exp[-- / d9A{9, -9) coth{9/ 2)] (81) 
2 J—c 



oo 



where A{9i,92) is defined in ([78|) . Since the exponents commute on constant, we can use 
the formula 

= (82) 

and obtain the following result: 
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< /i(X) >= exp( — / de{l + l/cosh^)e-2"^^°^'^^) 
8 JO 

= e^v{-\[Ko{2mX) + fC_i(2mX)]} (83) 
8 

< a{X) >= exp(-m|X|) exp{-l[Ko{2mX) + K_i{2mX)]} (84) 

8 

that is 

n{X) =< fi,{X)a2{X)as{X)fio{X) > 

= exp(-m|X|) exp{-l[3Ko{2mX) + Ko{6mX)3K_i{2mX) + K^i{6mX)]} (85) 

8 

This vacuum average behaves as X^^/^ at X << and decays exponentially at X >> 



VII. CONCLUSIONS AND DISCUSSION 



In this paper we have shown that doped spin-1/2 ladder systems are described by the 
random mass Majorana (Dirac) fermion model. On the basis of this model, we have calcu- 
lated the thermodynamic functions for the spin ladders. In particular, we predict 1/Tln^T 
low-temperature asymptotics for the linear magnetic susceptibility. This behavior is quite 
different from the simple Curie law. As discussed in Section |VD| , there is already good 
evidence of this behavior in numerical simulations. Of course, we hope that more precise ex- 
perimental measurements would be able to distinguish these nontrivial disorder effects from 
the contributions of uncorrelated 1/2 spins induced by impurities. We would like also to 
point out that the recent neutron data have shown the existence of the gap states, while 
the magnitude of the gap itself does not change with doping. This is certainly consistent 
with our theoretical results. 

We did not attempt to discuss more complicated questions related to the behavior of the 
correlation functions in such disordered systems. We only note here that the divergency of 
the density of states (and of the localization length) in the middle of the gap makes these 
systems different from standard one-dimensional disordered systems giving rise to a non- 



trivial critical regime at low energies p7| , |3T| . In fact, in the recent months, there has been 
quite an impressive progress in the understanding of the correlation functions. For instance, 
some important insight into the zero energy properties of the random mass Dirac model 
was provided by mapping of this problem onto a Liouville quantum mechanics [|l6l. An 



interplay between the critical regime at low energies and the standard localization regime 
was explored by Beresinskii's diagram technique Jl^ and by the supersymmetry method 

However, the influence of the disorder on the staggered magnetic susceptibility in spin 
ladders has been poorly studied so far. This quantity is important from the experimental 
point of view. It is vital for the understanding of the apparent promotion of the antiferro- 
magnetic ordering upon doping, which was experimentally observed in both the spin ladder 
and the spin-Peiels systems [^p8|. It is our opinion that future theories of the antiferromag- 
netic transition in these systems will be based on the mapping onto the random mass fermion 
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model, presented in this paper, in conjunction with the theoretical progress in dealing with 
such fermionic models achieved in Refs HTBI-ITsl . 
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FIGURES 



The QMC simulation data of the doping effect on uniform magnetic susceptibihty in 
spin ladders of Ref. []1T| are replotted as a function of l/(lnT)^ in comparison with Eg. (p9|) . 
The temperature T is measured in units of J. The fitting formula is xT = c [ a + 
6/(lnT)2], where a = 0.185(3), 0.145(2), 0.126(1), and b = 0.43(3), 0.29(1), 0.23(1) for 
doping concentrations c = 0.01, 0.05, 0.1, respectively. The dotted lines are expectations 
for uncorrelated free 1/2 spins induced by impurities, while the arrows on the left side 
indicate the renormalized values, anticipated from the renormalization group analysis ||13|. 
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